;********************************************
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;load library for  wetbulb_stull
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/heat_stress.ncl"
;load "./panel_two_sets.ncl" ; have to download this script to run
;********************************************
 begin

;set baseline
 dirPath = "~/source_data/"
 dirPath2 = dirPath+"/outdoor/"

 mo  = (/"ACCESS-CM2","AWI-CM-1-1-MR","BCC-CSM2-MR","CMCC-CM2-SR5","EC-Earth3","GFDL-ESM4", \
            "IITM-ESM", "KIOST-ESM","MIROC6","MIROC-ES2L","MPI-ESM1-2-HR","MRI-ESM2-0"/);

 case = (/"ssp585" /)
 ncase = dimsizes(case);
 nmod = dimsizes(mo)
 levs = (/1,1.5,2,2.5,3,3.5,4/); 7 warming levels
 nlev = dimsizes(levs)

 syr = 1990
 fyr = 2099;
 csyr = sprinti("%0.4i",syr);
 cfyr = sprinti("%0.4i",fyr);
 nyr = fyr - syr + 1;
 ndays = nyr*365

 ;read dimension for 1 deg data
  dimPath=dirPath+"/domain/"
  landfrac = "sftlf_fx_360x180_1deg_gn.nc"
  area = "areacella_fx_360x180_1deg_gn.nc"

  fdim  = addfile(dimPath+landfrac, "r")
  fdim2 = addfile(dimPath+area, "r")
  lat := fdim->lat
  lon := fdim->lon
  nlat = dimsizes(lat);
  nlon = dimsizes(lon);
  garea = fdim2->areacella  ;m^2
  lfrac = fdim->sftlf/100; %
  wgt = garea*1e-6*lfrac ; sqrkm
  wgt@_FillValue = lfrac@_FillValue; make sure FillValue is set correctly for wgt
  copy_VarCoords(garea, wgt)
  printVarSummary(wgt)
  print("global total land area: "+sum(wgt)*1e-6+"million sq km")



;===Use precalculated global warming amount relative to 1980-2009 (based on ERA5, HadCRUT5, BerkeleyEarth, NOAAGlobalTemp) to classify warming levels for each CMIP6 model
  ;see Data_global-warming-amount-ERA5-30yr-CMIP6_12models.ncl
  filename1 = dirPath+"Global_warming_amount_10yrma_"+syr+fyr+"_12models_obs30yr.csv"; 10-year running mean
  filename2 = dirPath+"Global_warming_amount_30yrma_"+syr+fyr+"_12models_obs30yr.csv"; 20-year running mean

;---Read in file as array of strings so we can parse each line
  lines  = asciiread(filename1,-1,"string")
  lines2 = asciiread(filename2,-1,"string")
  k = 1 ;number of header lines
  nlines = dimsizes(lines)-k   ; skip header lines
  delim = ","
;---Read fields (columns)
  models  =          str_get_field(lines(k:),1,delim)
  cases   =          str_get_field(lines(k:),2,delim)
  central_10yr =  tointeger(str_get_field(lines(k:),3,delim))
  period10yr =       str_get_field(lines(k:),5,delim)
 ;read mean temperature anomaly from file 1
  Ta_anomaly10yr = tofloat(str_get_field(lines(k:),7,delim))
  Ta_anomaly10yr@_FillValue=1e+20

  central_30yr =    str_get_field(lines2(k:),3,delim)
  period30yr =       str_get_field(lines2(k:),5,delim)
  Ta_anomaly30yr = tofloat(str_get_field(lines2(k:),7,delim))
  Ta_anomaly30yr@_FillValue=1e+20


 ;------choose whether 10-yr or 30-yr running period
  Ta_anomaly = Ta_anomaly30yr; either 10-yr or 30-yr running median/max
  periods = period30yr
  central_year = central_30yr

  levels := fspan(1,4,7); 1 to 4 degC by 0.5 degC
  nlev = dimsizes(levels)


  ;===mask G by specific subpopulations===
  ;--*read population distribution data from "Fig.S2.Population_map_water_rural.ncl"
  ;1) population spending >30mins to access water (rescaled to UNICEF,WHO gross values)
  ;2) rural farmer population data made from 1km SMOD and 1km POP data from EU GHSL
  ;as described in Supplementary Method 2 "Population distribution data"
  dirNc =  dirPath
  filNc2  = "population_data_watercollector_farmer.nc"
  pthNc2  = dirNc+filNc2
  ncdf2 = addfile(pthNc2 ,"r")
  POPR_1deg = ncdf2->POPR_1deg; population engaged in farming
  POPR_dens = ncdf2->POPR_dens
  POPW_1deg = ncdf2->POPW_1deg; population engaged in water collection
  POPW_dens = ncdf2->POPW_dens

  ;--mask global land area by water collection data
  water_area = where(POPW_1deg .gt. 0, wgt, wgt@_FillValue);
  print("Total area of water collection in 2020: "+sum(water_area)*100*1e-6 +" Mha")

  ;--mask global land area by rural population data
  rural_area = where(POPR_1deg .gt. 0, wgt, wgt@_FillValue); rural area
  print("Rural land area in 2020: "+sum(rural_area)*100*1e-6 +" Mha")



 PLOT  = True; 
 if (PLOT) then

   filNc = "G_3h_diel_localtime_outdoor.nc"
   pthNc  = dirPath2+filNc
   ncdf = addfile(pthNc ,"r")        ; open output netCDF file

   G_sun_med  = ncdf->G_sun_med
   G_sun_low  = ncdf->G_sun_low
   G_sun_high  = ncdf->G_sun_high
   G_diff_med  = ncdf->G_diff_med
   G_diff_low  = ncdf->G_diff_low
   G_diff_high  = ncdf->G_diff_high
   G_zero_med  = ncdf->G_zero_med
   G_zero_low  = ncdf->G_zero_low
   G_zero_high  = ncdf->G_zero_high
   ;levs = ncdf->lev
   Gday_sun_med  = ncdf->Gday_sun_med; [nlev,nlat,nlon]


   ;**use a common set of grid cells (those where Gday>0 emerges at 1C warming under the sun scenario) to calculate the diel profiles for all scenarios so that they are comparable
   ;-select grid cells with Gday>0 at 1C warming under sun scenario (~176 cells)
   Gday_4d := conform_dims(dimsizes(G_sun_med),Gday_sun_med(0,:,:),(/2,3/)) ; conform to 4d: [nlev,nstep,nlat,mlon] 

   G_sun_med := where(Gday_4d .gt. 0, G_sun_med, G_sun_med@_FillValue)
   G_sun_low := where(Gday_4d .gt. 0, G_sun_low, G_sun_low@_FillValue)
   G_sun_high := where(Gday_4d .gt. 0, G_sun_high, G_sun_high@_FillValue)
   G_diff_med := where(Gday_4d .gt. 0, G_diff_med, G_diff_med@_FillValue)
   G_diff_low := where(Gday_4d .gt. 0, G_diff_low, G_diff_low@_FillValue)
   G_diff_high := where(Gday_4d .gt. 0, G_diff_high, G_diff_high@_FillValue)
   G_zero_med  := where(Gday_4d .gt. 0, G_zero_med, G_zero_med@_FillValue)
   G_zero_low  := where(Gday_4d .gt. 0, G_zero_low, G_zero_low@_FillValue)
   G_zero_high := where(Gday_4d .gt. 0, G_zero_high, G_zero_high@_FillValue)


   ;---further mask by areas of subpopulation (water collection or farming)
   Amask_4d := conform_dims(dimsizes(G_sun_med),water_area,(/2,3/)) ; conform to 4d: [nlev,nstep,nlat,mlon] 
   ;Amask_4d := conform_dims(dimsizes(G_sun_med),rural_area,(/2,3/)) ; conform to 4d: [nlev,nstep,nlat,mlon] 
   G_sun_med := where(Amask_4d .gt. 0, G_sun_med, G_sun_med@_FillValue)
   G_sun_low := where(Amask_4d .gt. 0, G_sun_low, G_sun_low@_FillValue)
   G_sun_high := where(Amask_4d .gt. 0, G_sun_high, G_sun_high@_FillValue)
   G_diff_med := where(Amask_4d .gt. 0, G_diff_med, G_diff_med@_FillValue)
   G_diff_low := where(Amask_4d .gt. 0, G_diff_low, G_diff_low@_FillValue)
   G_diff_high := where(Amask_4d .gt. 0, G_diff_high, G_diff_high@_FillValue)
   G_zero_med  := where(Amask_4d .gt. 0, G_zero_med, G_zero_med@_FillValue)
   G_zero_low  := where(Amask_4d .gt. 0, G_zero_low, G_zero_low@_FillValue)
   G_zero_high := where(Amask_4d .gt. 0, G_zero_high, G_zero_high@_FillValue)

  
   print("")
   count_sun = num(.not.ismissing(G_sun_med(0,0,:,:))); count number of true values
   print("number of grid cells with Gday(sun)>0 at 1C: "+count_sun+" (masked by water/rural population)")
   print("")

   ;--calculate mean diurnal pattern of 3-hourly G across grid cells where annual maximum Gday>0
   ;NOTE: the 8 time steps are uniform GMT time, corresponding to different time zones in different grid cells
   G0_sun_diurnal_med = dim_avg_n(G_sun_med, (/2,3/)) ; [nlev,nstep]
   G0_sun_diurnal_low = dim_avg_n(G_sun_low, (/2,3/)) 
   G0_sun_diurnal_high = dim_avg_n(G_sun_high, (/2,3/)) 
   G0_diff_diurnal_med = dim_avg_n(G_diff_med, (/2,3/)) ; [nlev,nstep]
   G0_diff_diurnal_low = dim_avg_n(G_diff_low, (/2,3/)) 
   G0_diff_diurnal_high = dim_avg_n(G_diff_high, (/2,3/)) 
   G0_zero_diurnal_med = dim_avg_n(G_zero_med, (/2,3/)) ; [nlev,nstep]
   G0_zero_diurnal_low = dim_avg_n(G_zero_low, (/2,3/)) 
   G0_zero_diurnal_high = dim_avg_n(G_zero_high, (/2,3/)) 

   ;--select the warming amount 
   levid := get1Dindex(levels, 1)
   G0_sun_diurnal_med1 = G0_sun_diurnal_med(levid,:)
   G0_sun_diurnal_low1 = G0_sun_diurnal_low(levid,:)
   G0_sun_diurnal_high1 = G0_sun_diurnal_high(levid,:)
   G0_diff_diurnal_med1 = G0_diff_diurnal_med(levid,:)
   G0_diff_diurnal_low1 = G0_diff_diurnal_low(levid,:)
   G0_diff_diurnal_high1 = G0_diff_diurnal_high(levid,:)
   G0_zero_diurnal_med1 = G0_zero_diurnal_med(levid,:)
   G0_zero_diurnal_low1 = G0_zero_diurnal_low(levid,:)
   G0_zero_diurnal_high1 = G0_zero_diurnal_high(levid,:)
   levid := get1Dindex(levels, 2)
   G0_sun_diurnal_med2 = G0_sun_diurnal_med(levid,:)
   G0_sun_diurnal_low2 = G0_sun_diurnal_low(levid,:)
   G0_sun_diurnal_high2 = G0_sun_diurnal_high(levid,:)
   G0_diff_diurnal_med2 = G0_diff_diurnal_med(levid,:)
   G0_diff_diurnal_low2 = G0_diff_diurnal_low(levid,:)
   G0_diff_diurnal_high2 = G0_diff_diurnal_high(levid,:)
   G0_zero_diurnal_med2 = G0_zero_diurnal_med(levid,:)
   G0_zero_diurnal_low2 = G0_zero_diurnal_low(levid,:)
   G0_zero_diurnal_high2 = G0_zero_diurnal_high(levid,:)
   levid := get1Dindex(levels, 3)
   G0_sun_diurnal_med3 = G0_sun_diurnal_med(levid,:)
   G0_sun_diurnal_low3 = G0_sun_diurnal_low(levid,:)
   G0_sun_diurnal_high3 = G0_sun_diurnal_high(levid,:)
   G0_diff_diurnal_med3 = G0_diff_diurnal_med(levid,:)
   G0_diff_diurnal_low3 = G0_diff_diurnal_low(levid,:)
   G0_diff_diurnal_high3 = G0_diff_diurnal_high(levid,:)
   G0_zero_diurnal_med3 = G0_zero_diurnal_med(levid,:)
   G0_zero_diurnal_low3 = G0_zero_diurnal_low(levid,:)
   G0_zero_diurnal_high3 = G0_zero_diurnal_high(levid,:)
   levid := get1Dindex(levels, 4)
   G0_sun_diurnal_med4 = G0_sun_diurnal_med(levid,:)
   G0_sun_diurnal_low4 = G0_sun_diurnal_low(levid,:)
   G0_sun_diurnal_high4 = G0_sun_diurnal_high(levid,:)
   G0_diff_diurnal_med4 = G0_diff_diurnal_med(levid,:)
   G0_diff_diurnal_low4 = G0_diff_diurnal_low(levid,:)
   G0_diff_diurnal_high4 = G0_diff_diurnal_high(levid,:)
   G0_zero_diurnal_med4 = G0_zero_diurnal_med(levid,:)
   G0_zero_diurnal_low4 = G0_zero_diurnal_low(levid,:)
   G0_zero_diurnal_high4 = G0_zero_diurnal_high(levid,:)
    
   opt = True
   opt@PrintStat = True
   ;cstat = stat_dispersion(G0_sun_diurnal_med, opt) ;detailed statistics
   ;set color range between lower/upper 5% and min/max
   print("Min-max range of G0_sun_diurnal: ")
   printMinMax(G0_sun_diurnal_high,0)

   print("G0_sun_diurnal at 4degC: "+G0_sun_diurnal_med4(:))
   print("G0_sun_diurnal_low at 1degC: "+G0_sun_diurnal_low1(:))
   print("G0_sun_diurnal_high at 1degC: "+G0_sun_diurnal_high1(:))


  ; draw barplot and 25-75th percentiles around the ensemble mean
   nstep = 8 ; timesteps per day
   xx = ispan(0,21,3); +1.5

   xp    = new( (/2*nstep/), float )
   yp_s1   = new( (/2*nstep/), float )
   yp_d1   = new( (/2*nstep/), float )
   yp_z1   = new( (/2*nstep/), float )
   yp_s2   = new( (/2*nstep/), float )
   yp_d2   = new( (/2*nstep/), float )
   yp_z2   = new( (/2*nstep/), float )
   yp_s3   = new( (/2*nstep/), float )
   yp_d3   = new( (/2*nstep/), float )
   yp_z3   = new( (/2*nstep/), float )
   yp_s4   = new( (/2*nstep/), float )
   yp_d4   = new( (/2*nstep/), float )
   yp_z4   = new( (/2*nstep/), float )

   xp(0:(nstep-1)) = xx
   xp(nstep:(2*nstep-1)) = xx(::-1) ; reverse

  ;ensemble 25-75th percentile
   yp_s1(0:(nstep-1)) = G0_sun_diurnal_high1 ; uppper bound
   yp_s1(nstep:(2*nstep-1)) = G0_sun_diurnal_low1(::-1) ; lower bound
   yp_d1(0:(nstep-1)) = G0_diff_diurnal_high1 ; uppper bound
   yp_d1(nstep:(2*nstep-1)) = G0_diff_diurnal_low1(::-1) ; lower bound
   yp_z1(0:(nstep-1)) = G0_zero_diurnal_high1 ; uppper bound
   yp_z1(nstep:(2*nstep-1)) = G0_zero_diurnal_low1(::-1) ; lower bound
   yp_s2(0:(nstep-1)) = G0_sun_diurnal_high2 ; uppper bound
   yp_s2(nstep:(2*nstep-1)) = G0_sun_diurnal_low2(::-1) ; lower bound
   yp_d2(0:(nstep-1)) = G0_diff_diurnal_high2 ; uppper bound
   yp_d2(nstep:(2*nstep-1)) = G0_diff_diurnal_low2(::-1) ; lower bound
   yp_z2(0:(nstep-1)) = G0_zero_diurnal_high2 ; uppper bound
   yp_z2(nstep:(2*nstep-1)) = G0_zero_diurnal_low2(::-1) ; lower bound
   yp_s3(0:(nstep-1)) = G0_sun_diurnal_high3 ; uppper bound
   yp_s3(nstep:(2*nstep-1)) = G0_sun_diurnal_low3(::-1) ; lower bound
   yp_d3(0:(nstep-1)) = G0_diff_diurnal_high3 ; uppper bound
   yp_d3(nstep:(2*nstep-1)) = G0_diff_diurnal_low3(::-1) ; lower bound
   yp_z3(0:(nstep-1)) = G0_zero_diurnal_high3 ; uppper bound
   yp_z3(nstep:(2*nstep-1)) = G0_zero_diurnal_low3(::-1) ; lower bound
   yp_s4(0:(nstep-1)) = G0_sun_diurnal_high4 ; uppper bound
   yp_s4(nstep:(2*nstep-1)) = G0_sun_diurnal_low4(::-1) ; lower bound
   yp_d4(0:(nstep-1)) = G0_diff_diurnal_high4 ; uppper bound
   yp_d4(nstep:(2*nstep-1)) = G0_diff_diurnal_low4(::-1) ; lower bound
   yp_z4(0:(nstep-1)) = G0_zero_diurnal_high4 ; uppper bound
   yp_z4(nstep:(2*nstep-1)) = G0_zero_diurnal_low4(::-1) ; lower bound


 imagePath = dirPath+"../image"
 image1 = imagePath+"/Fig.7"; focus on the subpopulation engaged in water collection
 ;image1 = imagePath+"/Fig.S11"; focus on the subpopulation engaged in farming

  wks1 = gsn_open_wks("pdf",image1)
 
  plot1 = new(3,graphic)
  plot2 = new(3,graphic)
  plot3 = new(3,graphic)
  plot4 = new(3,graphic)
  res = True
  res@gsnDraw              = False
  res@gsnFrame             = False
 
  res@vpXF                 = 0.21
  ;res@vpWidthF             = 0.6
  ;res@vpHeightF            = 0.35
  res@vpWidthF             = 0.34  ;plot width
  res@vpHeightF            = 0.3  ;plot length
 
  res@tiXAxisFontHeightF   = 0.014;0.022
  res@tiYAxisFontHeightF   = 0.014;0.022
  res@tmXBLabelFontHeightF = 0.012
  res@tmYLLabelFontHeightF = 0.012
  res@tmYRLabelFontHeightF = 0.012
  res@gsnStringFontHeightF = 0.014
  res@tmXTOn               = False
  res@tmXBMinorOn          = False
  res@tmYLMinorOn          = True
  res@tmYRMinorOn          = False
  res@tiYAxisJust          = "BottomCenter"
  ;res@tiYAxisFontAspectF   = 1.5
  res@pmLegendDisplayMode    = "Never";
 

  res@tmYROn              = False
  res@tmYRLabelsOn        = False
  res@trYAxisType = "LinearAxis"
  res@gsnStringFontHeightF = 0.014;
  res@gsnLeftStringOrthogonalPosF = 0.03 ;adjust the space/distance of gsnLeftString/gsnCenterString
  res@gsnCenterStringOrthogonalPosF = 0.03
  res@tiYAxisJust      = "CenterCenter"
  ;res@tiYAxisOffsetXF  = -0.01


 ;-----------plot1----------
 res@vpYF          = 0.90   ; location of plot1
 res@vpXF          = 0.10    ; location of plot1
  res@gsnLeftString       = "~F22~a";
  res@gsnCenterString  = "1 ~S~o~N~C warming"; "Grid cells with G~S~day~N~>0"
  res@tiYAxisString    = "G (W m~S~-2~N~)"
  res@tiXAxisString    = ""; "Time"
  res@trXMinF          = 0
  res@trXMaxF          = 24
  res@tmXBMode         = "Explicit"
  res@tmXBValues       = xx
  res@tmXBLabels       = xx

  res@trYMinF          = -520;
  res@trYMaxF          := 400;
  res@tmYLOn              = True
  res@tmYLMode            = "Manual"
  res@tmYLTickStartF      = -500
  res@tmYLTickEndF        = 400
  res@tmYLTickSpacingF    = 100
  res@tmYLMinorPerMajor   = 4
  res@tmYLAutoPrecision = False
  res@tmYLPrecision     = 3; decimal place

  ;=====G under zero solar radiation (dark)=========
  res@xyLineColors      = "slategray"; "deepskyblue2"
  res@xyDashPatterns    = 0; 3
  plot1(0) = gsn_csm_xy(wks1, xx, G0_zero_diurnal_med1, res); zero solar
  gsres                   = True                        ; poly res
  gsres@gsFillOpacityF    = 0.1
  gsres@gsFillColor     = "slategray4";                ; color chosen
  dummy13 = gsn_add_polygon (wks1,plot1(0),xp,yp_z1,gsres);

  res@tiXAxisString    = " "
  res@tiYAxisString    = " "
  res@gsnLeftString    = " ";
  res@gsnCenterString  = " ";
  res@tmXBOn              = False
  res@tmYLOn              = False
  res@tmYLLabelsOn        = False
  res@tmXBLabelsOn        = False

  ;=======G under full solar radiation (sun)===========
  res@xyLineColors          = "purple2";
  plot1(1) = gsn_csm_xy(wks1, xx, G0_sun_diurnal_med1, res);
  gsres@gsFillColor       = "purple";                ; color chosen
  dummy11 = gsn_add_polygon (wks1,plot1(1),xp,yp_s1,gsres) ;sun

  ;=======under only diffuse radiation (shade)=======
  res@xyLineColors        = "orange2"
  plot1(2) = gsn_csm_xy(wks1, xx, G0_diff_diurnal_med1, res);
  gsres@gsFillColor       = "orange"                    ; color chosen
  dummy12 = gsn_add_polygon (wks1,plot1(2),xp,yp_d1,gsres) ;shade

  ;=add 0 reference line
  yvalues = (/0,0/)
  res@gsnYRefLine      = yvalues
  res@gsnYRefLineColor = "grey"
  res@gsnYRefLineDashPattern = 1
  res@gsnYRefLineThickness   = 1
  res@pmLegendDisplayMode    = "Never"
  plot01 = gsn_csm_y(wks1,yvalues,res) ; use same res as base plot
  overlay(plot1(0),plot01); plot01 becomes part of plot1



 ;------------plot2---------
 ;use vpXF to place two plots horizontally
  res@vpXF          = 0.6-0.05       ; location of plot2
  res@gsnLeftString       = "~F22~b";
  res@gsnCenterString  = "2 ~S~o~N~C warming"; "Grid cells with G~S~day~N~>0"

  res@tiYAxisString    = "G (W m~S~-2~N~)"
  res@tiXAxisString    = ""; "Time"
  res@trXMinF          = 0
  res@trXMaxF          = 24
  res@tmXBMode         = "Explicit"
  res@tmXBValues       = xx
  res@tmXBLabels       = xx
  res@tmXBOn              = True
  res@tmYLOn              = True
  res@tmYLLabelsOn        = True
  res@tmXBLabelsOn        = True

  res@trYMinF          = -520;
  res@trYMaxF          := 400;
  res@tmYLOn              = True
  res@tmYLMode            = "Manual"
  res@tmYLTickStartF      = -500
  res@tmYLTickEndF        = 400
  res@tmYLTickSpacingF    = 100
  res@tmYLMinorPerMajor   = 4
  res@tmYLAutoPrecision = False
  res@tmYLPrecision     = 3; decimal place


  ;=====G under zero solar radiation (dark)=========
  res@xyLineColors      = "slategray"; "deepskyblue2"
  res@xyDashPatterns    = 0; 3
  plot2(0) = gsn_csm_xy(wks1, xx, G0_zero_diurnal_med2, res); dark
  gsres                   = True                        ; poly res
  gsres@gsFillOpacityF    = 0.1
  gsres@gsFillColor     = "slategray4";                ; color chosen
  dummy23 = gsn_add_polygon (wks1,plot2(0),xp,yp_z2,gsres);

  res@tiXAxisString    = " "
  res@tiYAxisString    = " "
  res@gsnLeftString    = " ";
  res@gsnCenterString  = " ";
  res@tmXBOn              = False
  res@tmYLOn              = False
  res@tmYLLabelsOn        = False
  res@tmXBLabelsOn        = False

  ;=======G under full solar radiation (sun)===========
  res@xyLineColors          = "purple2";
  plot2(1) = gsn_csm_xy(wks1, xx, G0_sun_diurnal_med2, res);
  gsres@gsFillColor       = "purple";                ; color chosen
  dummy21 = gsn_add_polygon (wks1,plot2(1),xp,yp_s2,gsres) ; sun

  ;=======under only diffuse radiation (shade)=======
  res@xyLineColors        = "orange2"
  plot2(2) = gsn_csm_xy(wks1, xx, G0_diff_diurnal_med2, res);
  gsres@gsFillColor       = "orange"                    ; color chosen
  dummy22 = gsn_add_polygon (wks1,plot2(2),xp,yp_d2,gsres) ; shade

  ;=add 0 reference line
  plot01 = gsn_csm_y(wks1,yvalues,res) ; use same res as base plot
  overlay(plot2(0),plot01); plot01 becomes part of plot2


 ;-----------plot3----------
 res@vpYF          = 0.50   ; location of plot3
 res@vpXF          = 0.10    ; location of plot3
  res@gsnLeftString       = "~F22~c";
  res@gsnCenterString  = "3 ~S~o~N~C warming"; "Grid cells with G~S~day~N~>0"

  res@tiYAxisString    = "G (W m~S~-2~N~)"
  res@tiXAxisString    = "Hour of day"
  res@trXMinF          = 0
  res@trXMaxF          = 24
  res@tmXBMode         = "Explicit"
  res@tmXBValues       = xx
  res@tmXBLabels       = xx
  res@tmXBOn              = True
  res@tmYLOn              = True
  res@tmYLLabelsOn        = True
  res@tmXBLabelsOn        = True

  res@trYMinF          = -520;
  res@trYMaxF          := 400;
  res@tmYLOn              = True
  res@tmYLMode            = "Manual"
  res@tmYLTickStartF      = -500
  res@tmYLTickEndF        = 400
  res@tmYLTickSpacingF    = 100
  res@tmYLMinorPerMajor   = 4
  res@tmYLAutoPrecision = False
  res@tmYLPrecision     = 3; decimal place

  ;=====G under zero solar radiation (dark)=========
  res@xyLineColors      = "slategray"; "deepskyblue2"
  res@xyDashPatterns    = 0; 3
  plot3(0) = gsn_csm_xy(wks1, xx, G0_zero_diurnal_med3, res); dark
  gsres                   = True                        ; poly res
  gsres@gsFillOpacityF    = 0.1
  gsres@gsFillColor     = "slategray4";                ; color chosen
  dummy33 = gsn_add_polygon (wks1,plot3(0),xp,yp_z3,gsres);

  res@tiXAxisString    = " "
  res@tiYAxisString    = " "
  res@gsnLeftString    = " ";
  res@gsnCenterString  = " ";
  res@tmXBOn              = False
  res@tmYLOn              = False
  res@tmYLLabelsOn        = False
  res@tmXBLabelsOn        = False

  ;=====G under full solar radiation (sun)===========
  res@xyLineColors          = "purple2";
  plot3(1) = gsn_csm_xy(wks1, xx, G0_sun_diurnal_med3, res);
  gsres@gsFillColor       = "purple";                ; color chosen
  dummy31 = gsn_add_polygon (wks1,plot3(1),xp,yp_s3,gsres) ; sun

  ;=======under only diffuse radiation (shade)=======
  res@xyLineColors        = "orange2"
  plot3(2) = gsn_csm_xy(wks1, xx, G0_diff_diurnal_med3, res);
  gsres@gsFillColor       = "orange"                    ; color chosen
  dummy32 = gsn_add_polygon (wks1,plot3(2),xp,yp_d3,gsres) ; shade

  ;=add 0 reference line
  plot01 = gsn_csm_y(wks1,yvalues,res) ; use same res as base plot
  overlay(plot3(0),plot01); plot01 becomes part of plot3


 ;------------plot4---------
 ;use vpXF to place two plots horizontally
  res@vpXF          = 0.6-0.05       
  res@gsnLeftString       = "~F22~d";
  res@gsnCenterString  = "4 ~S~o~N~C warming"; "Grid cells with G~S~day~N~>0"

  res@tiYAxisString    = "G (W m~S~-2~N~)"
  res@tiXAxisString    = "Hour of day"
  res@trXMinF          = 0
  res@trXMaxF          = 24
  res@tmXBMode         = "Explicit"
  res@tmXBValues       = xx
  res@tmXBLabels       = xx
  res@tmXBOn              = True
  res@tmYLOn              = True
  res@tmYLLabelsOn        = True
  res@tmXBLabelsOn        = True

  res@trYMinF          = -520;
  res@trYMaxF          := 400;
  res@tmYLOn              = True
  res@tmYLMode            = "Manual"
  res@tmYLTickStartF      = -500
  res@tmYLTickEndF        = 400
  res@tmYLTickSpacingF    = 100
  res@tmYLMinorPerMajor   = 4
  res@tmYLAutoPrecision = False
  res@tmYLPrecision     = 3; decimal place

  ;=====G under zero solar radiation (dark)=========
  res@xyLineColors      = "slategray"; "deepskyblue2"
  res@xyDashPatterns    = 0; 3
  plot4(0) = gsn_csm_xy(wks1, xx, G0_zero_diurnal_med4, res); dark
  gsres                   = True                        ; poly res
  gsres@gsFillOpacityF    = 0.1
  gsres@gsFillColor     = "slategray4";                ; color chosen
  dummy43 = gsn_add_polygon (wks1,plot4(0),xp,yp_z4,gsres);

  res@tiXAxisString    = " "
  res@tiYAxisString    = " "
  res@gsnLeftString    = " ";
  res@gsnCenterString  = " ";
  res@tmXBOn              = False
  res@tmYLOn              = False
  res@tmYLLabelsOn        = False
  res@tmXBLabelsOn        = False

  ;=====G under full solar radiation (sun)=========
  res@xyLineColors          = "purple2";
  plot4(1) = gsn_csm_xy(wks1, xx, G0_sun_diurnal_med4, res);
  gsres@gsFillColor       = "purple";                ; color chosen
  dummy41 = gsn_add_polygon (wks1,plot4(1),xp,yp_s4,gsres) ; sun

  ;=======under only diffuse radiation (shade)=======
  res@xyLineColors        = "orange2"
  plot4(2) = gsn_csm_xy(wks1, xx, G0_diff_diurnal_med4, res);
  gsres@gsFillColor       = "orange"                    ; color chosen
  dummy42 = gsn_add_polygon (wks1,plot4(2),xp,yp_d4,gsres) ; shade

  ;=add 0 reference line
  plot01 = gsn_csm_y(wks1,yvalues,res) ; use same res as base plot
  overlay(plot4(0),plot01); plot01 becomes part of plot4


  ;--add some text
  txres                     = True                 ; text mods desired
  txres@txFontHeightF       = 0.012                ; default size is HUGE!
  txres@txJust              = "CenterLeft"         ; puts text on top of bars
  ;gsn_text(wks1,plot1(0),"Hottest 1% of land grid cells", 1.2, res@trYMaxF*0.95,txres) ; on top

  ;--add legend for shade
  lbres                    = True
  lbres@lbBoxMajorExtentF  = 0.6           ; smaller value for more space between color boxes
  lbres@lbMonoFillPattern  = True
  lbres@lbPerimOn          = False
  lbres@lbBoxLinesOn       = False
  lbres@lbFillOpacityF     = 0.2; 0.1
  lbres@lbFillColors       = (/"slategray4","orange","purple"/)
  labels = (/"","",""/); no labels for the shade bars

  ;legend for lines, more general function gsn_legend_ndc
  lgres                    = True
  lgres@lgLineThicknessF   = 1.5                 ; line thicknesses
  lgres@lgPerimOn          = False               ; no perimeter
  lgres@lgLineColors       = (/"slategray","orange2","purple2"/)
  lgres@lgDashIndexes      = (/0,0,0/)          ; line patterns
  labels := (/"G~S~day~N~ (dark)","G~S~day~N~ (shade)","G~S~day~N~ (sun)"/)
  lgres@lgLabelOffsetF     = 0.08
  lgres@lgLabelFontHeightF = 0.07;          ; font height (also scales with vpWidthF)

  labels := (/"","",""/); no labels for the shade bars
  lbres@vpWidthF           = 0.18          ; labelbar width
  lbres@vpHeightF          = 0.08;0.10          ; labelbar height
  lbres@lbFillColors       := (/"slategray4","orange","purple"/)
  gsn_labelbar_ndc(wks1,3,labels,0.07,0.9,lbres) ; draw left labelbar column
  lgres@lgDashIndexes      := (/0,0,0/)          ; line patterns
  lgres@vpHeightF          = 0.08              ; height of legend (NDC)
  lgres@vpWidthF           = 0.11              ; width of legend (NDC)
  lgres@lgLabelFontHeightF = 0.07;          ; font height (also scales with vpWidthF)
  lgres@lgLineColors       := (/"slategray","orange2","purple2"/)
  ;labels = (/"G~S~day~N~ > 0 (dark)","G~S~day~N~ > 0 (shade)","G~S~max~N~ > 0 (sun)"/)
  labels := (/"Dark","Shade","Sun"/)
  gsn_legend_ndc(wks1,3,labels,0.12,0.9,lgres)

 draw(plot1)
 draw(plot2)
 draw(plot3)
 draw(plot4)
 frame(wks1)


 end if




end 
